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Abstract 

The method of alternation projections (MAP) is an iterative procedure for finding the pro- 
jection of a point on the intersection of closed subspaces of an Hilbert space. The convergence 
of this method is usually slow, and several methods for its acceleration have already been 
proposed. In this work, we consider a special MAP, namely Kaczmarz' method for solving 
systems of linear equations. The convergence of this method is discussed. After giving its ma- 
trix formulation and its projection properties, we consider several procedures for accelerating 
its convergence. They are based on sequence transformations whose kernels contain sequences 
of the same form as the sequence of vectors generated by Kaczmarz' method. Acceleration can 
be achieved either directly, that is without modifying the sequence obtained by the method, 
or by restarting it from the vector obtained by acceleration. Numerical examples show the 
effectiveness of both procedures. 

1 Introduction 

Let Qi denote the orthogonal projection on a closed subspace Mj of an Hilbert space H, and let 
Qm be the composition of the r projecting operators Qj, that is Qm = Qr ■ ■ -Qi- Let Pm be the 
projection on M, the intersection of the subspaces Mj. 

We are looking for the projection of a given point x on M. It holds 

lim Ql^^t = Pm^, Vx e H. 

n— ^oo 

The method of alternating projections (MAP) consists in the iterations 

x„+i = Qa/x„, n = 0, 1,..., Xo=x. 
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This method often converges quite slowly and it needs to be accelerated [22]. For this purpose, the 

projection operator Qm can be replaced by another (non necessarily linear) operator T which can 
depend on n. 

There are two ways of using T 

1. Keep the sequence (x„) as given by MAP unchanged, and consider the new sequence (y„) 
built by y„ = Tx„. 

2. Iterate the operator T, that is consider the new iterations x„+i = Tx„. 

There exist many choices for T (the norms are the Euclidean ones). For example, the choice 



with = (x„, x,„ - Qmx„)/||x,„ - Qmx„P was proposed in |22l p. 44], while t„ = (r„, v„)/ (v„, v„) 
where r„ = Qm^u — and v„ = Q\^^n — "^Qm^ti + was discussed in [27] . More choices are 
presented in [19] . 

Kaczmarz's method [M] for solving a system of linear equations was proposed on 1937. Later, 
it was rediscovered by Gordon et al. [29], and applied in medical imagining. They called it ART 
[Algebraic Reconstruction Technique) , and its original version, or some variants, continue to be used 
for tomographic imaging. It is a particular case of row projection methods which received much 
attention (see, for example, [3|Bt[TT | [T9 | B7] ) . and it enters into the framework of MAP. It is also well 
suited for parallel computations and large-scale problems, since each step only requires one row of 
the matrix A (or several rows simultaneously in its block version), and no matrix-vector products. 
For an impressive list of publications on Kaczmarz's method, see [18]. Kaczmarz's method is also 
often used for solving an overdetermined consistent M x N linear system with M > N, but, in 
this paper, we restrict ourselves to the case of a square regular system. It is easy to extend the 
algorithms and the theory by properly replacing by M. 

A drawback of Kaczmarz's method, as in general of all projection iterative methods, is its often 
slow convergence. Thanks to its matrix analysis, we will be able to show how its convergence can 
be accelerated by some particular choices of the operator T corresponding to the two ways of using 
it described above. 

For definitions and properties of projections, see, for example, [11]. In the sequel, the same 
notation will be used for a matrix and the projection it represents. 

2 Kaczmarz's method 

We consider a. N x N linear system Ax = b. One single step of Kaczmarz's method consists in 



Tx. 



-n 




Pn+1 — Pn 



+ 



(b - Apn,ei) ^^j, 



(1) 
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where ej is the ith vector of the canonical basis of . There exist several strategies for choosing 
the index i at each step. The most common one is i = n (mod. A^) + 1. In this case, the method 
is called the cyclic Kaczmarz's method or, simply, the Kaczmarz's method, since it was originally 
proposed by Kaczmarz under this form [53] . It corresponds to a restarting from the result obtained 
after N single steps, that is, in other words, to a renumbering of them or, again in other words, to 
the extraction of a subsequence. Thus, one iteration of Kaczmarz's method consists in a complete 
cycle of steps in their natural order, that is 



Po = 



(b - Api_i,ei) J, 



= PN 

Denote by = A^ei the column vector formed by the zth row of A. Thus, the computation of 
each vector pj in ([2]) does not require any matrix-vector product since (b — Api^i, ej) = (b, e^) — 
[pi-i, A^ei). Thus, if we denote by bi the zth component of the right hand side b, then the 
computation of Pi is simply given by 

Pi = Pi-i H u — [75 aj. (3) 

This remark is one of the main advantages of Kaczmarz' method, and it allows an easy parallel 
implementation. 

Let Mi = {y I (b — Ay, ej) = 0}. Then is the oblique projection of pj_i on Mi along A'^ei. 
Moreover, setting 

, _ (b - Api_i,ei) 



one iteration of Kaczmarz's method ([2]) writes 

x„+i = X, + A^Ar,, A„ = (Ai, . . . , XnV e R^. 

Remark 1 

Let Yn he the iterates obtained by applying the GaussSeidel method to the system AAl'-j = b. Then 
X = A^y and x„ = A^yn J^'- 

Let us now analyze each step of ([2]). Inside one iteration, we have the following orthogonality 
properties for i = 1, . . . , (see [20ll26]) 



(b - Api,ei) = 0, 

(pi - pi„i,x - Pi) = 0, 
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and it follows 



= ||x - pif + 2Ai(b - Api,ei) + ||pi - Pi_i|p, 

since (b — Api, ej) = 0. Thus 

iix - Pi-if = iix - Piii^ + iipi - Pi_l||^ 

which shows that ||x — pj|| < ||x — pj_i||. Therefore, ||x — x„_|_i|| < ||x — x„||. Thus, as proved 
in [2B], these inequalities are strict, and Kaczmarz's method is always converging to the solution 
of the system. 

Obviously (b — ylpj_i,ej)^ < (b — Api^i,h — Api^i). In the case where (b — Apj_i,ej)^ > 
(b — Api_i,b — Api-i)/N, one can prove, by an analysis similar to what is done in |20l p. 122], 
that the following result holds 

Proposition 1 

// (b - Api_i, Bif > (b - Api_i, b - Api_i)/iV, then 

I|x-P.f <(l-^^;^^) ||x-p._,f. 

Then, in the general case, this coefficient could be an approximation of the convergence factor 
of the method. This result is quite similar to a corresponding result proved in [32] for Gastinel's 
method |25] and other methods, or to an extension (Thm. 4.27 in [23]) of another result given 
in [32]. 

2.1 Matrix interpretation 

Following [26], where it seems that it first appeared, let us give the matrix interpretation of Kacz- 
marz's method (ED (see also [I|[20|H7]). 

We set 

Pi 

Q^ 

Qi 

We have 

Qi = A~'PiA = I 



= I - ActieJ, 
= A~'P,A, 
= h-Ap,. 

_ A^BiejA ^ a,af ^ 
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The matrix Pi represents the obhque projection on along AA^Si, while Qi is the rank — 1 or- 
thogonal projection on (A^Bi)^ along A^ei. Thus, for any vector y, (P^y, e^) = and {Qiy, A^ei) = 
0. 

Inside one iteration, we have, for i = 1, . . . ,N, 

Pi = QiPi-i + (b, ei)a.i = QiPi-i + (/ - Qi)x 
x-pi = x-pi„i - (£>i„i,ei)Q:i = (5i(x-pi„i) } (4) 
Qi = Qi-i - iOi-i, fii)AoLi = PiQi_i. 

The first relation can be written 

Pi = QiP^^l + A-\I - Pi)h. 

Notice that A~^{I - Pj)b = (b, ei)ai. 

The second relation in (j4]) can also be written as 

X - p, = X - Pi 1 - a^ef ^,„i, 

and, replacing oti by its expression, it follows 

II Il2 II ||2 c^iSi-l^^i) /■ aT 

||x-pi|| = ||x-pi_i|| - 2 (x- pi_i, A e 

12 (^i-l'S 



\2 



X - Pi_l 



2 

2 II 



Summing up this identity for i = 1, . . . , A^, we obtain an expression for the gain of one iteration 
of Kaczmarz's method 



N , 

I l|2 _ II ||2 \Qi-l^ 

|x — X^^^lll — ||X — X„|| — 



\2 



Setting 



i=l " *ll 



P = Pn---Pi 

Q = Qn---Qi= A~'Pn---P,A = A^'PA 

TTn = b AXn, 



it holds 

X - X„+i = Q(x - x„) 
The matrix Q represents an orthogonal projection, but not P. It obviously follows 



(5) 



x-x„ = g"(x-Xo), r„, = P"ro, n = 0, 1, . . . , (6) 
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and we have 



X„+i = Q^n + A~\I - P)h. 



(7) 



Since, after the first iteration of Kaczmarz's method (which ends with Pn), we set xi = p^r, and 
we start the second iteration from po = xi, it means that we are continuing the original steps ([1]), 
and that the new vector pi is, in fact, the vector Pat+i of ([1]), that the new vector p2 is PAr+2, and 
so on. Thus, Kaczmarz's method is only a renumbering of the single steps, as previously explained, 
keeping only those whose index is a multiple of A^, that is x„+i = p(„_|_i)7v and r„+i = Q(^n+i)N- The 
main interest of this renumbering lies in the relations ([5]) which express the connection between two 
consecutive iterations by means of the fixed matrix Q, instead of relations using matrices changing 
at each step of an iteration. We similarly have Q(^n+i)N+j~i = Q^^^ QnN+j-i^ 3 = l;---;-^! with 
Q'^'J^ = Pj-i ■ ■ ■ PiPn ■■■Pj- Notice that Q^^^ is identical to Q. 

The matrices P and Q are similar, and it holds p{Q) < 1 as proved in [26]. Thus, the sequence 
(x„) generated by Kaczmarz's method converges to the solution x of the system. Moreover, it 
follows from standard results on iterations of the form ([6]) that 

||x-x„|| = o(p(g)"). 

By slightly improving Thm. 4.4 of [23], which is based on a result by Meany [35] on the norm 
of a product of orthogonal projections of rank — 1, we have the 

Proposition 2 

iix-x ,ir^< /'i-_(^fi^'liix-x 11= 
" [ n.liM'-e.pj" 

A generalization of this result was recently given in [2] for the case where A and b are partitioned 
into blocks of rows. 

Let us give some details about the block version. Assume that the matrix A is partitioned into 
the blocks of rows Aj, A^, . . ., where the matrix Aj G M^'^^ contains the rows A''i + ■ ■ ■ + Ni^i + 1 
up to A''i + ■ • • + A'j.i + Ni of A (with A^^q = 0), and the vector b is partitioned accordingly. Each 
single step of the method now consists in the treatment of a block as a whole, and one complete 
cycle of all the blocks in their natural ordering is called the (cyclic) block Kaczmarz's method. As 
for Kaczmarz's method, we are able to give a matrix interpretation of this extension which remains 
valid with now cti = (AJY = Ai{AjAi)-^ and Pi = I - AonEj , where Ei G M^""^" is the matrix 
whose columns are the vectors eN^+...+Ni^i+i) ■ ■ ■ , eAr-^+..._|_7Vi_i+AriOf the canonical basis of M^, and 
it holds Qi = I — Ai{Aj Ai)~^ Aj . With these notations, one step of the block Kaczmarz's method 
writes 

p, = p,_i + (Af)tEf(b-Ap,_i). 
This relation clearly generalizes (|3]). 
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3 Convergence acceleration 



For accelerating the convergence of a sequence, it can be transformed into another one which, 
under some assumptions, converges faster to the same hmit. The idea behind such a sequence 
transformation is to assume that the sequence to be accelerated behaves as a model sequence, 
or satisfies some property, depending on unknowns parameters. The set of these sequences is 
called the kernel of the transformation. The unknown parameters are determined by imposing 
that the sequence interpolates the model sequence from a certain starting index n. Then, the 
limit of the model sequence, which depends on n since the parameters depend on it, is taken as 
an approximation of the limit of the sequence to be accelerated which is thus transformed into 
a new sequence. Of course, by construction, if the initial sequence belongs to the kernel of the 
transformation then, for all n, its limit is obtained. Although this observation was never proved 
rigourously, if the sequence is close in some sense to the kernel of the transformation, there is a 
good chance that it will be accelerated. This is why the notion of kernel is so important. 

In practice, when having to accelerate a given sequence, one can use a known sequence trans- 
formation (also called an acceleration algorithm or an extrapolation method), and verify that it can 
be accelerated by it. Another approach is, starting from some algebraic property of the sequence 
to be accelerated, to construct a special transformation adapted to it. In both cases, the behavior 
of the sequence has to be analyzed or has to be characterized by some property. 

On convergence acceleration methods for vector sequences, see, for instance, [T3] . 

3.1 The sequence generated by Kaczmarz's method 

Consider the sequence of vectors (x„) obtained by Kaczmarz's method. Let Hi, be the minimal 
polynomial of the matrix Q for the vector x — xq, that is the polynomial of smallest degree v < N 
such that n^(Q)(x - xq) = 0. Since n^(Q)(x - xq) = v4-in^(P)A(x - xq) = A-^T{^{P)ro, is 
also the minimal polynomial of P for the vector yq. Setting ny(^) = cq + ci^ + ■ ■ ■ + Cy^^ , it holds 
from ([6]) 

Q'TI,,((5)(x - xo) = co(x - x„) + ci(x - x„+i) H h c,,(x - x„+,,) = 0, n = 0, 1, . . . (8) 

Thus, it follows that the vectors x — x„ produced by Kaczmarz's method satisfy such an homoge- 
neous linear difference equation of order z/, whose solution is well-known. 

If Q is nondefective and if we denote by ri, . . . , the distinct zeros of 11;^, this solution writes 



where the cij's are scalars and the Vj's the eigenvectors corresponding to the r^'s. If Q is defective, 
the expression for x — x„ still involves powers of the eigenvalues, but it is more complicated. 




= 0, 1, . . . , 



i=l 
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3.2 Sequence transformations 



The vector x can be exactly computed if we are able to build a sequence transformation whose 
kernel consists of sequences of the form ([8]). However, since, in practical applications, v could be 
quite large, we will restrict ourselves to building a sequence transformation whose kernel contains 
all sequences of the form 



with k < u. Thus, solving this equation for the vector x, gives us an approximation of it denoted 
y^"^ since it depends on k and n. 

For that purpose, one has to find a procedure for computing the unknown coefficients a^^^ , • • • , a^^^ , 
taking into account that they are numbers while x — x„_|_j are vectors. All such transformations can 
be considered as vector generalizations of Shanks' transformation [lOj for sequences of numbers. 
The common idea behind these transformations is to obtain a system of linear equations whose 
solution is Cq"'' , . . . , a^"'' , and then to compute y["^ ~ x. 

It turns out that several vector sequence transformations based on (or including) such a kernel 
already exist and have been studied by various authors: the various e-algorithms [71I1H1I19], the 
Minimal Polynomial Extrapolation (MPE) [17], the Modified Minimal Polynomial Extrapolation 
(MMPE) [711371I15], the Reduced Rank Extrapolation (RRE) [211I36], Germain-Bonne transforma- 
tions [28], the if-algorithm [16], and the i?-algorithm [10]. They are described and analyzed, for 
example, in [5l[71[8l[l0lllllll2lllll6]. 

Writing ([9]) for the indexes n and n + 1, subtracting, and multiplying scalarly by a vector y 
leads to the scalar equation 



where A is the difference operator defined by Ax„ = x^+i — x„. 

As explained in [HI p. 39], there are several possible strategies for constructing our system 
which, apart from an additional normalization condition, has to contain k equations 

• use only one vector y, and write ffTOj) for the indexes . . . ,n + k — 1, 

• write (ITOl) only for the index n, and choose k linearly independent vectors y, 

• write several relations fllOp . and choose several linearly independent vectors y. 

Adding to these k equations the condition Cq"^ + ■ ■ ■ + a^"'' = 1, which does not restrict the 
generality (and is needed since the sum of the coefficients must be nonzero in order for x in (I9l) to 
be uniquely defined), we obtain a system of A; + 1 equations in the k + 1 unknowns ao"\ . . . , a^"^ 
(which depend on n and k). Its coefficients are denoted and, for any of the preceding strategies, 
this system writes 




x„) + ar(x - x„+i) H h 4 - ^n+k) = 0, = 0, 1, . . . , 



(9) 



ar(y, Ax„) + af\y, Ax„+i) + ■ ■ ■ + ^"^(y, Ax„+fc) = 0, 



(10) 




I = 



1 
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where the coefficients d!fj, for i 



k and j = 0, . . . , k, will be given below according to the 

(n 
k 



chosen strategy. Then, for a fixed value of k, our sequence transformation (x„) i — (yl"^) is defined 
by 



which can be written as 



(n) 



(n) 



(n) 



n 



0,1, 



(12) 



(n) 



"1,0 



"fc,0 



1 

(n) 



"i,o 



in) 
k,0 



^n+k 



1 

(n) 



(n) 



Xj^ AXn 

"i,o ""i,o 



d 



in) 
k,0 



6d 



in) 
k.O 



Ax„+fc_i 



""fc,A;-l 



(5d 



(n) 

l,fc-l 



(13) 



where 5 is the difference operator defined by Sd^^ 
expression shows that the vector y^""* is the Schur complement [11] 



d-l^j+i — dfj ■ The second determinantal 



(n) 

yl 



x„ [/^x„, . . . , Ax„ 



(n) 
fc,0 



541i / 



However, y[."^ is not a projection. 

The second determinantal formula in (|T3i) means that all our sequences transformations can be 
also written as 

y^") = x„ - «S"^Ax„ al")Ax„+,_i, = 0, 1, . . . , (14) 

where the a^^^^'s are solution of the linear system 



(15) 



Let us now specify various choices of the coefficients d^^j , for i = 1, . . . , k and j = 0, . . . , k 



3.2.1 The vector Shanks' transformations 

We first consider sequence transformations which are directly inspired by the scalar sequence trans- 
formation of Shanks [IQ], and can be recursively implemented by Wynn's e-algorithm [18] or its 
generalizations, or by the ^'-algorithm [10], or the if-algorithm [T6] . 
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Choosing d!^^ = (ep, Ax„_|_j+j_i), and replacing Xj in the first row of the numerator of the first 



ratio in ( 1131) by (ep,Xj) corresponds to Shanks' transformation applied componentwise to the 
vector sequence (x„), and it gives the pth component of the vector y["^ This transformation 
can be recursively implemented, separately for each component p = 1, . . . ,N, hj the scalar 
e-algorithm of Wynn 



• The e-algorithm can be also applied directly to a sequence of vectors after defining the inverse 
u^^ of a vector u as u^^ = u/(u, u). This vector e-algorithm was also proposed by Wynn |49j . 
However, for this algorithm, ( |T3l) is no longer valid and the determinants have to be replaced 
either by bigger and more complicated ones [30], or by designants which generalize them 
in a non-commutative algebra [38l[39]- Thus, in this case, an underlying system of linear 
equations for y^""* does not exist, and this transformation has to be recursively implemented 
by the vector e-algorithm whose rules are, for a sequence (u„) of real (to simplify) vectors in 

e^^l = OgM^, 4"^ = u„gM^, n = 0,l,... 

Jn) _ (n+1) / (n+1) {n)x_l u r, 

tfc+l — £k-l + l^fc ) ; K, n — U, i, . . . 

• The choice d\"j = (y, Ax„+j+j_i), where y is any nonzero vector so that the denominators 
of ( |T3l) differs from zero, leads to the topological Shanks' transformation introduced in [7]. 
It can be implemented via the topological e-algorithm. In this algorithm, the inverses are 
defined in a different way for an even or an odd lower index as 

(4r'^ - 4? = y/(y, 4r'^ - 4:^) 

/ (n+1) _ (n) X 1 _ / (n+1) _ (n)x , ( Jn+1) _ Jn) (n+1) _ Jn). 
l^2fe+l ^2k+l) — \^2k ^2kJ/\^2k+l ^2fc+l; ^2fe ^2k}- 

In the preceding e-algorithms, only the vectors (or the numbers) with an even lower index e^^J 
are interesting for the purpose of convergence acceleration, those with an odd lower index s^^_^-j^ 
being intermediate computations. Applying any of them to the sequence (sq"^ = x„) produced 
by Kaczmarz's method gives £^2k ~ yl""*- Notice that the computation of one vector e^^k requires 
the 2k + 1 vectors x„, . . . ,^n+2k- Thus, when k is increased by 1, the computation of each new 
vector needs two additional iterates of Kaczmarz's method, while it only requires one when n 
increases. 

3.2.2 The MPE, MMPE, and RRE 

The following choices also lead to known transformations but, in this case, the computation of y^""* 
only requires the k + 2 vectors x„, . . . , x^+^+i. Algorithms for their recursive implementation also 
exist IE 
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The choice ci-"-^ = [ji, Ax„+j), for z = 1, . . . , /c and j = 0, . . . ,k corresponds to the MMPE 
(Modified Minimal Polynomial Extrapolation) [71I37II45] , where the are linearly independent 
vectors. This transformation is another generalization of the topological Shanks transforma- 
tion also given in [7]. 



• If we take d[^j = (Ax„+i_i, Axn+j), for i = 1, . . . , A; and j = 0, . . . , fc, we obtain the MPE 
(Minimal Polynomial Extrapolation) [17]. This method is mathematically equivalent to a 
transformation due to Germain-Bonne [28] . Other choices are proposed in the same reference. 

• The choice d\^j = (A^x„+j_i, Ax„+j) leads to the RRE (Reduced Rank Extrapolation) [2T|[36] . 

(n) 

It must be noticed that the vector always exists for the RRE but not for the MPE or the 
MMPE. Existence conditions are discussed in [33] and 1431. 



3.2.3 The simplest transformations 

For k = 1, the transformations have the following very simple forms 

• the MMPE and the topological Shanks' transformation write 

(„)_ (yi.AxQ 
" " (yi,A2xO 

• for the MPE and the transformation due to Germain-Bonne, we have 

(n) (Ax„, Ax„) 



(Ax„,A2xO 



for the RRE, it holds 



(n) _ _ (A^x„,Ax„) ^ 

yi — x„ Ax^, 

(A^x„, A^x„) 
which is the same as the one proposed in [27] . 

for the vector e-algorithm, we obtain 

(n+l) _ (n) . 



y) = x„ . 1 H — — —-T with £] 

J- "^-^ / (n+l) (n) (n+l) (n)\ J- 



_ - si"^) ' (Ax„, Ax„) 
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3.2.4 Orthogonality properties 

Using the determinantal formula (fT3|) . we have the following orthogonality properties 
Proposition 3 

The vector Q(y["'* — x) — (y^"-* — x) is orthogonal to 

• Ax„,, . . . , Ax„_|_fc_i for the MPE (which is identical to Germain-Bonne transformation) , 

• y for the topological Shanks transformation (topological e -algorithm), 

• A^x„, . . . , A^x„+fc_i for the RRE, 

• yi,...,yfc for the MMPE. 

For a discussion of the properties of these algorithms and their application to the solution of 
systems of linear and nonlinear equations, see [33] . 

4 Acceleration of Kaczmarz's method 

Although they are different, all the acceleration methods presented in Section |3] have the same 
kernel, namely sequences satisfying (j9]). Thus, when applied to any sequence of vectors generated 
by iterations of the form x„_|_i = -Bx„ + c, n = 0, 1, . . ., we obtain yi"^ = (J — B)~^c for all n, where 
V is the degree of the minimal polynomial of B for the vector x — xq. This result means that these 
methods are direct methods for solving systems of linear equations [7]. Therefore, this property 
holds for the iterates of Kaczmarz's method, since 11,^(1) ^ (otherwise the matrix I — Q would be 
singular). Of course, in practice, these algorithms cannot be used for obtaining the exact solution 
when the dimension of the system is large since, for computing yl"'', they require the storage of too 
many vectors. 

4.1 The AK and RK algorithms 

The preceding algorithms can be used in two different ways for accelerating the iterations (|2]) 
produced by Kaczmarz's method as it will now be explained. Since the computation of y[."'' does not 
require the same number of vectors issued from Kaczmarz's method, according to the transformation 
used, in order to have a unified presentation, we will denote by l+l the number of vectors Xj needed 
to compute y^""* by any of the preceding procedures. Remember that i depends on k and that i = 2k 
for the e-algorithms, and i = k + 1 for the other algorithms. 

• The first one consists to apply one of the algorithms implementing a sequence transformation 
to the sequence xq, xi, . . . given by Kaczmarz's method, and, after fixing the index k, to build 
simultaneously the sequence zq = y^^^ zi = y^^'', . . .. The computation of y^^'* can only begin 
after having computed the iterate x^, but the computation of each new transformed vector 
needs only one new iterate of Kaczmarz's method. This procedure is called the accelerated 
Kaczmarz (AK) algorithm. Let us give the general structure for the implementation of the 
AK algorithm. 
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Accelerated Kaczmarz (AK) algorithm 



Require A G M^^^, b G M^, xq G M" 
Choose A; G N, A; > 1 
Set i = k + loi (i = 2k 
for n = 0, 1, . . . until convergence do 
Po ^ x„ 

Compute Pi, i = 1, . . . , 

x„+i ^ Pat 

If > £ - 1 then 

Compute y^*^^^^^ 

(n-l+l) 

end if 
end for n 



• In the second way, we set Xq, and we compute Xi, . . . , x^ by Kaczmarz's method, we apply one 
of the algorithms implementing a sequence transformation to them, and we obtain zq = y^*^^ . 
Then, we restart Kaczmarz's method from Zo, that is we set xq = zq = y^^^ we compute the 
new ^ vectors xi, . . . , x^ by Kaczmarz's method, we apply again the acceleration algorithm 
to these vectors xq, . . . , x^, we obtain zi = yf \ we restart Kaczmarz's method from xq = zi, 
and so on. This second procedure is called the restarted Kaczmarz (RK) algorithm and it can 
be implemented as follows. 



Restarted Kaczmarz (RK) algorithm 



Require A G R^^^, b G M^, xq G M" 
Choose A; G N, /c > 1 
Set £ = A: + 1 or £ = 2A; 
for n = 0, 1, . . . until convergence do 
for J = 0, 1 do 

Po ^ 

Compute Pi, i = 1, . . . ,N 
Xj+i ^ Pat 
end for j 

Compute y^*^^ 
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Xo z„ 
end for n 



4.2 The case of the £— algorithms 

Let us now consider the particular case of the e-algorithms. 

If we assume that the eigenvalues of the matrix Q, defined in Section 12. ![ are numbered such 
that |ri| > |r2| > ■ ■ ■ > |r7v| > 0, then the sequences {s^2k)^ ioi k < u fixed, constructed by these 
algorithms are such that [6|H2|H6] 

||x-s£)||=0(r,V), n^oo. 

If the preceding assumptions on the Tj's are not satisfied (which corresponds to several eigenvalues 
of Q having the same modulus or to a defective matrix Q), the expression for x — x„ is more 
complicated than given above [H] (see [12] for the complete expression in the scalar case), but the 
convergence is still accelerated by the ^-algorithms. These results also hold for the topological e- 
algorithm, independently of the choice of the arbitrary vector y occurring in this algorithm. Thus, 
each sequence {s^2k) obtained by the accelerated Kaczmarz algorithm converges to x faster than 
the preceding sequence (£2fc-2) when n tends to infinity. Quite similar results hold for the RRE 
and the MPE 03]. 

If the topological 5-algorithm is applied to a sequence (x„) obtained by any iterative method of 
the form x„_|_i = i?x„ + c, with the choice y = fq = b — Axq, then the vectors Sg^^ /c = 0, 1, . . ., are 
identical to those obtained by Lanczos' method (that is, for example, by the biconjugate algorithm 
or the conjugate gradient algorithm in the symmetric positive definite case) [H Thms. 4.1 and 4.2, p. 
186-7]. This property comes out from the determinantal expressions of the vectors produced by the 
topological 5-algorithm and by Lanczos' method. From ([7]), we immediately see that this is the case 
of the vectors given by Kaczmarz's method which correspond to B = Q and c = A~^{I—P)h. Thus, 
if the topological e-algorithm is applied, with y = fq, to the sequence (x„) given by Kaczmarz's 
method, then the sequence {e'^^) is exactly the sequence produced by Lanczos' method. Using it as 
explained in the restarted Kaczmarz algorithm corresponds to restarting Kaczmarz's method with 
the vector obtained by the Lanczos' method after k of its iterates. At each restart, the vector y 
has to be taken equal to the corresponding residual. 

5 Implementation 

The implementation of our convergence acceleration procedures can be realized by three different 
ways. The first one consists, for a fixed value of k, in solving the linear system ( ITT]) and using f fT2]) . 
The second way is to employ (ITT]) and the system (fT5]) . The third possibility is to apply, when 
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it exists, a recursive algorithm. Since, for the vector e-algorithm, no underlying system of linear 
equations is known, its implementation can only be realized through its recursive rules, given in 
Section I3.2.1I 

An important practical problem is the choice of k. On one side, the effectiveness of our pro- 
cedures seems to increase with k but, on the other side, the number of vectors to be stored also 
increases with it. Thus, if the system is large, the value of k has to be kept quite small since too 
many vectors will have to be stored. The numerical stability of the procedures has also to be taken 
into consideration when k increases. 

Since Kaczmarz's method often converges slowly, the quantities Ax„+j are small, and the el- 
ements of the linear systems to be solved for obtaining the coefficients of the linear combination 
giving y^"'* approach zero. Thus, whatever the method used, the linear systems ( ITTj) or f|T5|) are 
ill-conditioned even for small values of k, and rounding errors can degrade the results. We tried 
several ways for computing their solution but the influence on the vectors y^"''* was small. 

Instead of solving the systems, it is possible to use a recursive algorithm for the computation of 
the vectors y^""*. Usually, sequence transformations are more efficiently implemented via a recursive 
algorithm as those mentioned above, and, in our case, this kind of approach seems to give results 
less sensitive to rounding errors for the topological and the vector e-algorithms. 

The quantities computed by all the recursive algorithms for implementing our transformations 
can be displayed in a triangular array (or a lozenge one for the e-algorithms) , and y^,""* corresponds 
to the lowest rightmost element of the array. Its computation imposes to compute first all the 
preceding elements in the array. Thus, it could be much too expensive to store all these vectors if 
k is not quite small and if is large. Hopefully, it is possible to proceed in the array by ascending 
diagonals (that is, starting from yg""'''^'* = ^n+k, to compute y["^'^~"^'', . . . ,y^k^), and, thus, storing 
only one ascending diagonal. For more details on such a technique, see |13] . 

Another problem is the choice of the arbitrary vector y in the topological Shanks transformation, 
and of the vectors yi, . . . , y^ appearing in the MMPE (see [33] for a discussion about this choice). 
This is an unsolved problem and further studies will be necessary. 

In the literature, the system Ax = b to be solved by Kaczmarz' method is often replaced by 
the system Z^Ax = Dh, where D = diag(l/||aj||). Thus, each row vector of the new matrix DA 
will have a norm equal to 1. This preconditioning, of course, simplify the computation of Pi in ([3]), 
and we will use it. 

As a last comment, notice that, although some transformations are equivalent from the theo- 
retical point of view (for instance the MPE and Germain-Bonne transformation, or the MMPE 
and the if-algorithm) , they do not produce exactly the same numerical results. The same is true 
when we compare the results of a transformation implemented in the different ways described in 
this paper. 
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6 Numerical results 



Our numerical experiments were performed using Matlab® 7.11 and matrices coming out from the 
gallery set. The solution was set to x = (1, ... , 1)^ and the right hand side b was computed 
accordingly. We took random vectors with components uniformly distributed in [— 1,+1] for the 
arbitrary vector y in the topological e-algorithm, and for the vectors yj in the MMPE. 

Some figures show the Euclidean norms of the errors and, in order to compare the acceleration 
brought by each procedure, we also give the ratios of the norms of the errors between the iterate 
z„ obtained by the AK or the RK algorithm and the iterate of Kaczmarz' method with the highest 
index used in its construction (for AK), or the iterate with the highest index which would have 
been used if we have led the method continue without restarting it (for RK) , that is 

J'^""""",, (AK) and „ ~ ^" (RK). 

||x„+£ — X|| ||X(„+i)(^+l) — x|| 

We consider the parter matrix A, N = 1000, k,{A) ^ 4.2306, a Toeplitz matrix with singular 
values near vr. In the Figures [1] and [21 we compare Kaczmarz' method with the MMPE, MPE, 
RRE and the topological e-algorithm, implemented by solving the system (fTT]) . and the vector 
e-algorithm, respectively for the AK and RK algorithms, with k = 5. 



AK errors: parter, N=1 000, k=5 ^ AK ratios: parter, N=1 000, k=5 




Figure 1: AK algorithm: errors and ratios for parter matrix, = 1000, k = 5. 



In these figures we see that all methods achieve a good precision with an advantage for the 
vector e-algorithm. Moreover, its convergence is smoother. The ratios grow up because all methods 
almost stagnate when a good precision is attained while the error of Kaczmarz' method continues 
to decrease slowly. In particular, for RK algorithm, the vector e-algorithm attains its full precision 
after 4 iterations while the AK algorithm needs more iterations. For this example, the dominant 
eigenvalue of A is 0.8732178, and the second one is 0.3170877. Thus, according to the theoretical 
results of Section 14. 2[ a good acceleration is observed with k = 1 for all procedures. 
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RK errors: parter, N=1000, k=5 



RK ratios: parter, N=1000, k=5 



Kac 

- MMPE 
■ ' ' MPE 

- RRE 

- EPS-TOP 

— EPS-VEC 




■ MMPE 
MPE 
RRE 

■ EPS-TOP 

■ EPS-VEC 



Figure 2: RK algorithm: errors and ratios for parter matrix, = 1000, k = 5. 



The Figure E] shows the ratios for the clement matrix, k,{A) ~ 1.13145 x 10^^, a tridiagonal 
matrix with zero diagonal entries, again with = 1000 and k = 5. For this example, the MPE 
and the RRE coincide. Again the vector e-algorithm is the best. 

For toeppen, a pentadiagonal Toeplitz matrix, with A^ = 1000, A; = 5, we have the Figure HI 
Notice that the MMPE and the topological e-algorithm do not work well and that this behavior 
could be due to the choice of the vectors y and y^. In fact, these choices can affect the results 
in a quite serious way. For instance, these methods sometimes exhibit better convergence and 
acceleration with y = (1,...,1)^ and y^ = e^. With the RK algorithm, considering the first 
50 iterations for k = 8, the vector e-algorithm attains a ratio of 10~^ at iteration 20 and after 24 
iterations a division by zero occurs. The ratios for the RRE and the MPE have a minimum of 10~^° 
at iteration 25. The topological e-algorithm diverges from the beginning. The MMPE exhibits an 
erratic convergence and the ratio goes down to 10"'^ at the iteration 42, and then increases. 

We also tried our procedures on a bigger matrix. The results for lesp, a tridiagonal matrix 
with real, sensitive eigenvalues, A^ = 10000, k = 5, k,{A) ^ 6.9553 x 10^, with the AK algorithm 
are given in Figure |5l We see that Kaczmarz' method and the acceleration procedures all attain 
full precision after 90 iterations. Thus all ratios will grow up. However, the vector e-algorithm 
attains an error of less than 10~^^ after about 20 iterations while Kaczmarz's method has only an 
error of order 10"^ at the same iteration. 

An important point in any iterative method is to have a quite reliable stopping criterion. Usually 
such iterations are stopped by using the residuals. However such a computation will need a matrix- 
vector product and, in our case, one of the interests of Kaczmarz' method will be lost. Thus, since 
the results given by our acceleration procedures often stagnate when some precision is attained while 
those of Kaczmarz' continue to decrease, the iterations can be stopped as soon as the following 
ratios grow up significantly 
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AK ratios: clement, N=1000, k=5 




- MMPE 
■ ' ' MPE 

- ' RRE 

- EPS-TOP 

— EPS-VEC 




RK ratios: clement, N=1000, W=5 
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■ MIVIPE 
MPE 
RRE 

■ EPS-TOP 

■ EPS-VEC 



Figure 3: AK and RK algorithms: ratios for clement matrix, = 1000, k = 5. 



AK ratios: toeppen, N=1000, k=5 




- MIVIPE 
■ ' ' MPE 

- ' RRE 

- EPS-TOP 

— EPS-VEC 



RK ratios: toeppen, N=1000, k=5 



■ MMPE 
MPE 
RRE 

■ EPS-TOP 

■ EPS-VEC 




10 15 20 25 30 



Figure 4: AK and RK algorithms: ratios for toeppen matrix, A^ = 1000, k = 5. 



J' (AK) and „^ " „ (RK). 

An example with the vector e-algorithm is given in the Figure [6] (left: AK algorithm, right: RK 
algorithm). The solid line corresponds to the norm of the error and the dashed one to the ratio for 
the stopping criterion. 

However, it must be noticed that, in the case of the RK algorithm, this stopping criterion 
involves iterates of Kaczmarz' method that have not been computed and used in the acceleration 
procedures. Thus, it is not usable in practice. Since the quantities ||Az„|| usually decrease rapidly, 
the iterates can be stopped when it is small enough and begins to stagnate or even to grow up. 

In our examples, we also add a white noise between 10~^ and 10^^ to the vector b. The norm 
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AK errors: lesp, N=10000, k=5 




AK ratios: lesp, N=10000. k=5 



20 40 60 80 100 




■ MMPE 
MPE 
RRE 

■ EPS-TOP 

■ EPS-VEC 



20 40 60 80 100 



Figure 5: AK algorithm: errors and ratios for lesp matrix, = 10000, k = 5. 

AK error and stopping ratios: iesp, N=2000, k=3 RK error and stopping ratios: iesp, N=2000, k=3 





Figure 6: AK (left) and RK (right) algorithms: errors (solid) and stopping (dashed) ratios for lesp 
matrix, A^ = 2000, A; = 3. 



of the error of the results obtained by our acceleration procedures attains the level of the noise in 
most cases. 

For the matrix baart of the Matlab Regularization toolbox [31j of dimension 120 whose condi- 
tion number is 2.28705 x 10^^, an error of the form e = (5||b||u/A/iV where 6 is between 10~^ and 
10~^ and u is a vector whose components are random variables from a normal distribution with 
mean and standard deviation 1, was added to b. The norm of the error achieved with the vector 
e-algorithm goes down to 10"° ''^ for k = 1,2, and 3, which is a little bit better than the results 
obtained in [T5] . 

Let us mention that the stopping criterion given above only works correctly for small noises. 
Maybe, it because the vector b is not involved in it. 
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7 Conclusions 



From our numerical results, it seems that the vector e-algorithm is the best procedure for acceler- 
ating Kaczmarz' method. However, the recursive implementation of the other procedures has to be 
tested numerically to see if it leads to better and more stable results. In our numerical examples, 
we tried several values of k. Although, when the dominant eigenvalue of A is well separated, k = 1 
leads to quite a good acceleration, it seems that higher values produce better results in general. 
Anyway, the choice of k and of the vectors y and yj are important points which need to be studied 
more deeply. Other recursive algorithms, such as those developed by Germain-Bonne [2H], or the 
VTT and the BVTT jT3], not considered in this work, have also to be tried. 

The acceleration of the Symmetric Kaczmarz' and the Randomized Kaczmarz' methods, which 
are also sequential row-action methods that update the solution using one row of A at each step, of 
other methods of the MAP class, of the SIRT (Simultaneous Iterative Reconstruction Techniques) 
methods, has to be considered. Finally, applications to tomography and, in general, to image 
reconstruction have to be considered. 
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